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The quality of simulated hypersonic stagnation region heating on tetrahedral meshes is 
investigated by using a three-dimensional, upwind reconstruction algorithm for the inviscid 
flux vector. Two test problems are investigated: hypersonic flow over a three-dimensional 
cylinder with special attention to the uniformity of the solution in the spanwise direction 
and hypersonic flow over a three-dimensional sphere. The tetrahedral cells used in the 
simulation are derived from a structured grid where cell faces are bisected across the 
diagonal resulting in a consistent pattern of diagonals running in a biased direction across 
the otherwise symmetric domain. This grid is known to accentuate problems in both 
shock capturing and stagnation region heating encountered with conventional, quasi-one- 
dimensional inviscid flux reconstruction algorithms. Therefore the test problem provides 
a sensitive test for algorithmic effects on heating. This investigation is believed to be 
unique in its focus on three-dimensional, rotated upwind schemes for the simulation of 
hypersonic heating on tetrahedral grids. This study attempts to fill the void left by the 
inability of conventional (quasi-one-dimensional) approaches to accurately simulate heating 
in a tetrahedral grid system. Results show significant improvement in spanwise uniformity 
of heating with some penalty of ringing at the captured shock. Issues with accuracy near 
the peak shear location are identified and require further study. 


I. Nomenclature 


Bold face, lowercase variable names refer to vectors. Bold face, uppercase variable names refer to matrices. 


Roman symbols 


A 

E 

fx> 

H 


MW 


n X ' 

nAi,hA2,nA3 


P 

q 

R 

T 


area 

total energy 

flux in x' direction, [pU x >,puU x > + pn x , x ' , pvU x > + pn y , x > , pwU x ' + pn ZjX >, pU x 'H] T 
total enthalpy 

unit vectors in Cartesian coordinate directions ( x,y,z ), respectively 
molecular weight 

unit vector in x' direction, n XtX 'i + n VtX 'j + n ZtX 'k 
unit vectors orthogonal to faces Al, A 2, A3 
pressure 

conserved variables [p, pu, pv, pw, pE] T 
characteristic variable in direction x ' , dq x ' = H X 'dq 
eigenvector matrix for flux Jacobian, = R^/A^/R^/ 
temperature 
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Roman symbols, continued 

U, V, w 

u x , 

x',y',z' 

Greek symbols 


Cartesian velocity components 
velocity component in x' direction, un XtX 
principal directions across element 


vn v 


■ wn z 


a 

P 

A 

P 

O 


coefficient relative to face for defining derivative in an element, Eq. 30 
coefficient relative to node for defining derivative in an element, Eq. 8 
diagonal matrix of eigenvalues for flux Jacobian, = Rjj, 1 A^/R^/ 
density 

element volume 


Subscripts 

1,2,3 

1-2, 1-3, 2-3 
A 
k 
L 

LSq 

n 

R 


node indices in triangular element 

edge identification as function of node endpoints 

element id, centroid of element 

face oriented index defined by index of opposite node in the element 
left virtual node in Option 1 

computed by Least Squares of gradients in surrounding elements 

node oriented index in element 

right virtual node in Option 1 

function of direction x' , usually involving U x ' 


II. Introduction 

The simulation of hypersonic flows with fully unstructured (tetrahedral) grids has severe problems with 
respect to the prediction of stagnation region heating. 1,2 The problems arise for two reasons. 

First, good shock capturing in the hypersonic regime requires alignment of the grid with the captured 
shock. Alignment here means that all surfaces of a control volume in contact with the shock are either 
parallel to the discontinuity or orthogonal to the discontinuity. Any skewness present in a conventional 
upwind scheme (in which first-order reconstruction is a function of only two nodes on opposite sides of a 
shared face) results in a smearing of the flow which manifests itself as a non-physical distribution of entropy 
from streamline to streamline crossing the shock. In the stagnation region, there is very little dissipation of 
these entropy gradients as they pass from the shock to the boundary-layer edge. Consequently, the surface 
heating is very sensitive to these non-physical gradients. Structured grid systems are generally easy to align 
with the captured bow shock. NASA’s main computational aerothermodynamic simulation codes LAURA 3,4 
and DPLR 5 both take great pains to align structured grid with the captured bow shock to improve solution 
quality but accept the limitation that internal shocks (e.g. wake recompression, control surface compression) 
are computed with no special grid consideration. (In fact, sharp double cone code validation test problems 6 
revealed that extremely fine structured grids were required to achieve a grid converged solution in a problem 
where structured grid could not be easily aligned with internally reflected shocks and shear layer over a 
separation bubble.'’ 8 ) Tetrahedral cell topologies by their very nature will always have at least one surface 
that is skewed to the captured shock. 

Second, tetrahedral grids are not well suited to the preservation of symmetry and the biased tetrahedral 
grids used in the tests herein are particularly ill-suited for this purpose. Asymmetric elements tend to induce 
non-physical cross flows. Again, this problem is particularly evident in the stagnation region where the 
physical flow velocities themselves are small, the non-physical, grid-induced velocities are of corresponding 
magnitude, and the surface heating is a sensitive indicator of these problems. The combination of random 
jaggedness in the shock capturing process as well as non-physical cross flow velocities cross-couple to produce 
unacceptable surface heating predictions. 

Problems with heating can be overcome if one uses semi-structured grid (e.g. prisms) across the boundary 
layer and adapts the grid to the shock. Nompelis et. al. 9 show excellent heating results for the cylinder 
test problem 10 on families of grids where a prismatic grid was adapted to the shock and acceptable heating 
results when an unbiased (random face orientation), tetrahedral grid was adapted to the shock. In this 
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case, random jaggedness may disguise some issues with shock capturing on tetrahedral elements. Prismatic 
elements are relatively easy to generate on blunt body geometries and their superior performance with respect 
to aeroheating predictions have led to their use as standard practice in unstructured simulations. There is 
no better alternative with todays algorithms than to use prismatic elements to capture the boundary layer 
and bow shock. Note however that application of shock fitting in the context of unstructured grids may offer 
improved solution quality by bypassing the difficulties associated with capturing a strong bow shock. 11 

If one is doing a hypersonic simulation without any free shear layers or internal shocks the specialized 
application of prismatic elements at the body and bow shock is perfectly acceptable. If such flow structures 
exist - as they do in almost all interesting problems - then the accuracy of any algorithm must be questioned 
when it ignores special topological grid requirements where the features are harder to resolve. 

A challenge problem was identified 2, 10 in which a highly biased, tetrahedral grid is applied to the sim- 
ulation of a hypersonic flow over a 3D cylinder at VJ*, = 5000m/s, poo = 0.001kg/m 3 and T\ ' x = 200K in 
air (M^ ~ 17)- (In fact, other hypersonic free stream conditions are acceptable — the challenge here is to 
recover the proper spanwise uniformity within 1%). A successful simulation requires a constant spanwise 
simulation of heating, pressure, and shear and symmetric simulation of these quantities on the right and left 
sides of the cylinder. The grid bias accentuates any issues with the algorithm. To date there have been no 
successful simulations of this challenge problem on the given grid using conventional upwind formulations. 
A recent thesis on a related unstructured grid without bias (random orientation of tetrahedra in spanwise 
direction) produced excellent heating symmetry using a high-order Discontinuous Galerkin finite element 
method with a PDE based artificial viscosity. 12 

A flux reconstruction algorithm that is not constrained by the local orientation of the grid is thought 
to have greater potential for a successful resolution of the challenge problem. The finite element method 
noted above 12 has an inherently three-dimensional stencil in the evaluation of flux. Other approaches using 
more conventional, upwind, finite- volume formulations are known as “rotated upwind schemes” . These were 
originally generated to address problems in capturing shocks on structured grids that were oblique to the 
grid. 13,14 The underlying ideas for this class of algorithm have been adapted for use in three-dimensional 
reconstruction using the node-based unstructured solver FUN3D and are applied to the challenge problem. 

III. Multi-Dimensional Reconstruction Algorithm 

A. Control Volume 

A triangular grid system is presented in Fig. 1. The two-dimensional reconstruction algorithm is developed 
relative to this figure. The extension of the algorithm to three dimensions is presented in Appendix A. 

Assume a node based scheme as illustrated in Fig. 1 with ( x,y ) coordinates of all nodes (1,2,3, ...) 
given and dependent variables q at these same nodes to be determined by solution of conservation laws. 
The conservation laws are discretized relative to the dual control volume, represented by dashed lines in 
Fig. 1. The dual volume around node 3 passes through the centroids of all elements surrounding the node 
(A, B 1 C, D, E, F ) and the midpoints of all edges separating these elements. 

B. Reconstruction Overview - Principal (Rotated) Direction 

Consider a conventional, quasi-one-dimensional, first-order reconstruction of inviscicl flux across face AcB 
separating the dual volumes around nodes 1 and 3. This flux is defined as a function of information at nodes 
1 and 3 of edge c: 

f AcB = \ [f 3 + fi - Rc|A c |R“ 1 (q 3 - qi)] (1) 

Inviscicl fluxes are then defined using a loop over edges and the reconstruction direction is fully determined 
by the grid. 

The basic idea in multi-dimensional reconstruction is that: (1) reconstruction directions should be based 
on local flow characteristics and not be constrained by the grid; and (2) the flux components for any face will 
utilize a stencil employing all nodes of a surrounding element. Whereas node based schemes loop over edges 
when using quasi-one-dimensional reconstruction the multi-dimensional reconstruction focuses on a loop over 
elements. Consider the two-dimensional element defined by nodes (1,2,3) with centroid A in Fig. 1. The 
computation of inviscicl flux in this element employs the same infrastructure for computing viscous flux. 
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Figure 1. Schematic of node-based, unstructured grid system showing element centroids (uppercase let- 
ters), edge midpoints (lowercase letters), nodes (numerals), the dual volume (dashed lines), and the principal 
directions (dotted lines with arrowhead). 


A principal direction in the element, x', is defined by the direction of Vp, computed using the same 
Green Gauss evaluation of the derivatives as applied for the viscous terms. If density is locally constant 
the principal direction aligns with the local velocity computed as the average of nodal values. Inviscicl flux 
within the element will be computed along and orthogonal to the principal flow direction. 

The flux components are defined using nodal weights based on the Green-Gauss formulation for derivatives 
in the component directions. Thus, 


diA 
dx 1 


E fkAkUk,. 

A k— faces 


(2) 


where ffc is the average value of f across all nodes defining surface k of element A. With reference to Fig. 1 
this derivative can be written 


df ‘x\A 
dx 1 


OL 3,x'{^2,x' + fl,x') + a l,x'(f3,x' + ^2,x') + + f3,x') 


(3) 


where 


tlfc,x' 


Aknk, x ’ 

2Sl A 


and k is a face index identified by the number of the opposite node. The sum Ylk=faces a k,x' 
closed element. It is convenient to reorder and scale the coefficients as follows. 


dfx\A 

dx ' 


(ct2,x' A C«3,x')fl + ( a 3,x' + Ot l,a/)f2 + (cti,x' + OL 2,x')^3 
[/3l,x'fl + 02,x'f2 A /33, x 'f3] 


(4) 

0 for any 

(5) 

(6) 


Note that Ax' can be interpreted as a distance across the element in x' direction. It is defined 


1 

Ax' 


^ ' (l^n— l,x' A O n _(_i i x / |) 

n=nodes 


(7) 


where a cyclic indexing is assumed. The reordering and scaling yield the following relations for /3 n> x'- 


Pn,x' 

— Ax (tin— l,x / A O n ^-i j x / ) 

(8) 

^ ^ Pn,x' 
n 

= 0 

(9) 

E IAbx'1 

= 1 

(10) 
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Two options for computing i x / are developed. Both mirror the baseline quasi-one-dimension reconstruction 
using Roe’s averaging across nodes. The first option averages dependent variables at two virtual nodes 
used in the reconstruction. The second option utilizes actual nodes and applies a weighted average to the 
computed flux on surrounding edges. 


C. Option 1 - Virtual Node Averaging 

A right and left virtual state for computing f x > are computed as follows. 

ttR,a;' = 2 /* ' (|/3n,x'| + 0n,x')qn 

n=nodes 

Q.L,x' = 2 ^ ^ d/^n,®' | — Pn,x')Qn 


(11) 

(12) 


n=nodes 


fx ' ~ 2 i 


+ l'( < l)fl,x' Rx' \^x'\(dq x ' dqx',«m)J (13) 

where the left and right states of f are functions of the left and right states of q, respectively. The matrices 
Rj;/ and Aj,/ are computed as functions of the Roe’s average of the right and left states. 


dqx' 


(qR,x' qL,x') 

Rx' ^hn— nodes Pn,x ,( ln 


(14) 


dq R: x' — Rx' ^ ) (|/?ro,x'| T 0n,x')^ q.n,LSq 

n=nodes 

dc{L,x' ~ Ax ^ ^ ftn,x , ')^'Qn,LSq 


n=nodes 


dq x ',um = minmod 


2dq fl , x / ,2dq x > ,2dq L ^ , ~{dq R ^ + dq LtX >) 


(15) 

(16) 

(17) 


D. Option 2 - Weighted Average of Edges to Principal Node 

Select the largest of \0 n ,x'\ to identify the principal node for construction of f x > and express it as a function 
of the remaining coefficients. Assume node 3 is the principal node. Then 


dfx' = 0\,x'f\ + 02,x'f2 + 03, x' /*3 

= 01,x'fl + /?2,x'/2 — {01, x' + 02,x')fd, (18) 

= 01, x' {fl — f 3 ) + 02, x' {f 2 ~ / 3 ) 

The reconstructed flux in the direction x' , f x i, is computed as a weighted average of surrounding edges. 
Furthermore, it is noted that if any of the surrounding edges are parallel to x' then the weight of that edge 
will equal 1 and the weight of the other edges will equal zero. 

fx' = \01,x' |fl-3,x' + \02,x'\^2-3,x’ (19) 


where 


fl— 3 ,x' — 2 fl,#' 

+ f 3}X ' - sign(/3i jX /)R 1 7 3iX / |Ax_ 3i3 ;/ |(dqi_ 3iX / - 

- dqi- 3tX i ,, 

lim) 

(20) 

^ 2 — 3 ,x' = ^ ^ 2 ’ x ' 

+ f 3 } x' - sign(/J 2 ,x')R7-3,x'l A 2-3,x'|( rf q2-3,x' " 

- dq 2 - 3 , x \. 

lim) 

(21) 

dQl- 3 ,x' = 

Ri-3,x'(qi - q3) 



(22) 

d^ 2 - 3 ,x' = 

R2-3,x'( c l2 — R3) 



(23) 

dc{l — 3 ,x',lim = 

minmod 

2^q 1 ,a; / , 2 dqi-s iX ' , 2dq 3 ,*' , - (dqi ?x / 

+ dq 3 , x /) 


(24) 

dc{ 2 — 3 ,x' ,lim = 

minmod 

2 dc[ 2 ,x' •> 2 dc\ 2 — 3 ,x' i 2 dc[ 3 x > , — {dc[ 2 ,x' 

+ dq3.x') 


(25) 

dQn,x' — 

Ar'R / 

^ Qn,LSqH-x' 



(26) 
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E. Component Assembly on Faces 


Both option 1 (Eq. 13) and option 2 (Eq. 19) present algorithms to define the flux, f x /, in the principal 
direction x' . The same logic is applied to construct f y >. Within each two-dimensional, triangular element 
there are three partial faces. The partial faces Ac and Ab can be combined into an equivalent face A3 
separating node 3 from the other nodes in the element. In like manner Ac and Aa form A1 and Aa and 
Ab form A 2. The flux across each equivalent face is constructed from the component flux in the principal 
direction and orthogonal to the principal direction. Thus, for element with centroid A as shown in Fig. 1 
the flux through each equivalent face is 


f An — ^An,xAx' T ^An^yAy' (27) 

where An refers to face separating node n from other nodes in the element. 

F. Shock Capturing 

Both options 1 and 2 are observed to admit temperature undershoots in capturing strong bow shocks. This 
undershoot occurs for strong shocks even in the case of a first-order reconstruction where <iq, x -',; iTO = 0. 
Though every component of the reconstructed flux uses a baseline Total Variation Diminishing (TVD) 
algorithm the net formulation is not TVD. A floor on temperature is specified, but this practice causes the 
convergence to stall as temperatures are reset. Modifying the viscosity within the shock layer to produce a 
cell Reynolds number equal to 8 has been observed to help this problem at the expense of a thicker captured 
shock. The algorithm is still being tested. Preliminary results indicate improvement in solution quality 
across the bow shock but degradation of solution quality in vicinity of shock / boundary-layer interactions. 

G. Accuracy 

Tests to evaluate order of accuracy for heating and shear on sequentially refined tetrahedral grids have been 
impeded by an inability to converge to a sufficiently low error norm. Issues with shock capturing noted 
above are believed to be the primary cause of the stalled convergence. A related issue is a ringing in the 
computation of the principal direction - a slight change in the computed density gradient direction effects 
the reconstruction. This ringing is easily suppressed by freezing the principal direction when the solution is 
converged below a user specified norm. 

The concern here is that the midpoint of the virtual nodes does not in general coincide with the centroid 
of the element; consequently, the reconstructed flux is offset from the dual volume surface by some small 
distance within the element. The blunt body solutions that follow match the pressure and shock shape of 
second-order structured grid benchmarks. However, there are differences in the shear and heating relative to 
structured grid benchmarks that are larger than expected for a second-order algorithm. Algorithms to create 
a second-order correction to the reconstructed flux that places it exactly at the centroid have not yielded 
any significant improvement. In fact, these “corrections” destroy the symmetry that is the focus of the 
test problems which follow. Thus, the current algorithm should be considered a work in progress. It offers 
significant improvement over the baseline quasi-one-dimensional reconstruction algorithm for hypersonic flow 
simulation on tetrahedral grids but there are still accuracy issues in the simulation of shear and heating. 

IV. Numerical Results 


A. FUN3D 

The three-dimensional reconstruction scheme is tested within FUN3D, which is a node based, fully unstruc- 
tured, finite volume solver of the Euler and Navier-Stokes equations. 15 The method uses Least Squares 
(LS) gradient information to execute second-order accurate inviscid flux reconstruction. 16 Viscous gradients 
are computed from a Green-Gauss formulation. A suite of modules includes all of the gas physics mod- 
els in LAURA and VULCAN 17 for thermodynamics, transport properties, chemical kinetics, and thermal 
relaxation. The baseline inviscid flux reconstruction algorithms utilize quasi-one-dimensional (edge-based) 
reconstruction using Roe’s averaging with Yee’s symmetric total variation diminishing algorithm adapted 
for unstructured grids. 
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B. Challenge Problem - Biased, Tetrahedral Grid on 3D Cylinder 

Accurate simulation of stagnation region heating in hypersonic flows is a key requirement for acceptance of 
any algorithm proposed for use in aerothermodynamic analyses. A structured grid solution generated with 
LAURA is used as both a benchmark and to generate initial grids for use in FUN3D as described originally 
in an earlier paper . 2 The test problem 10,18 uses Vx = 5000m/s, px = 0.001kg/m 3 , and Tx = 200K. 
Sutherland’s law for air is used to define transport properties in all perfect-gas cases. As noted previously, 
this simple problem provides insight into the ability of a scheme to cleanly capture the bow shock, smoothly 
resolve the post-shock stagnation region flow and predict a smoothly varying heating distribution around the 
stagnation point. These flowfield characteristics are particularly sensitive to the inviscid flux reconstruction 
algorithm and problems that are not evident in well aligned, structured (hexahedral) grids are exposed in 
the unstructured (tetrahedral) environment . 18 

The structured grid, adapted to the shock and boundary layer, from LAURA is converted from hexahedral 
elements to tetrahedral elements by adding diagonal edges consistently from minimum index to maximum 
index corners. A comparison of the grids in a plane of nodes perpendicular to the cylinder axis is presented 
in Fig. 2. The structured grid has 65 nodes from the body to the inflow boundary ahead of the bow shock 
and 61 nodes from the left to right outflow boundaries. Node placement is identical between the two grids. 
The placement of additional edges in the unstructured grid is the only difference. The strong biasing of 
diagonals in the grid is an intentional characteristic to expose algorithm deficiencies that may otherwise be 
averaged out in the simulation. 



Figure 2. Structured grid (LAURA) and biased, unstructured grid (FUN3D) in plane orthogonal to cylinder 
surface. 


A key element of this test problem is the addition of ten spanwise cells, shown in Fig. 3, across the 
cylinder, providing additional degrees of freedom in the simulation to allow asymmetries to develop. Earlier 
tests 18 (when the code was referred to as High Energy Flow Solver Synthesis - HEFSS) have shown that the 
single spanwise cell grids show good agreement with the structured code results in heating. (Fig. 4). The 
spanwise degrees of freedom enable non-physical cross-flow velocities to develop and irregular shock capture 
in the spanwise direction that combine to corrupt the predicted heating distribution. 

1. Structured Grid Reference Results 

Results from a previous paper 2 are briefly repeated here to show the symmetry qualities being sought in the 
unstructured formulation. Structured, hexahedral grids are assessed in FUN3D to confirm that spanwise 
and circumferential symmetry is preserved. A perfect gas case (7 = 5/3, MW = 2 M x = 4.25) is tested. 
A constant spanwise value of heating and shear is recovered when eigenvalue limiting is engaged as shown 
in Fig. 5. The near perfect overplotting of symbols in Fig. 5(a) for each of 11 spanwise nodes (ten spanwise 
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(a) Structured, quadrilateral surface elements 



(b) Unstructured, triangular surface elements 


Figure 3. Surface mesh on cylinder with ten rows of spanwise cells. 



Figure 4. Surface heating over cylinder at standard test conditions with 5-species reacting air model and 
fully catalytic wall. 


cells) at each x location around the cylinder indicates that the baseline algorithm is behaving as expected 
for the structured grid case and provide a metric by which to evaluate the unstructured results. 

2. Conventional Upwind Results on Biased, Tetrahedral Grid 

The perfect gas test case described above was repeated on the equivalent unstructured grid with identical 
location of nodes. The structured and unstructured results for heating and shear were compared in an earlier 
paper 2 where an approximate spanwise variation of ± 10 % about the mean was observed in peak values of 
heating and shear. There was also approximately 10% difference in the peak values of shear on the right 
and left sides of the cylinder for the unstructured grid. The large value of 7 = 5/3 and low molecular 
weight MW = 2 produce a relatively large sound speed and moderate supersonic Mach number Mra = 4.25. 
Even for this relatively low supersonic Mach number heating and shear asymmetries were produced on the 
challenge grid. 

Selecting a perfect gas with 7 = 7/5 and MW = 28.8 yields a lower sound speed and higher free stream 
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(a) Heating (blue squares) and shear (red squares) distribu- 
tion 



Figure 5. Results on cylinder computed on structured grid for 7 = 5/3 and MW = 2 with ten spanwise cells. 


Mach number Mqq = 17.34. This more severe condition produces even larger spanwise variation in heating 
and shear (Fig. 6(a)) with ±30% about the mean - an unacceptable result for aerothermodynamic analyses. 
Recall that the symbols at every x location should overplot and the peak value of shear on the right should 
equal the peak value of shear on the left. The LAURA result in this case is 52 W/cm 2 which is in good 
agreement with the spanwise mean of the heating at the stagnation point. The surface heating contours 
in Fig. 6(b) indicate that a slight drift in the velocity - probably induced by grid bias - produces a higher 
heating toward the front {y = 0) plane. 



(a) Heating (blue squares) and shear (red squares) distribu- 
tion for each of ten rows of spanwise cells. 



Figure 6. Simulation results for 7 = 7/6 and MW = 28.8 using conventional reconstruction on the challenge 
grid with ten spanwise cells. 


3. Option 1 Results on Biased, Tetrahedral Grid 

Option 1 results for the 7 = 7/5 and MW = 28.8 conditions as in the previous section are shown in Fig. 7. 
The Option 1 results are presented as symbols and a benchmark, grid-converged, structured grid solution 
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(LAURA) for the same case are presented as solid lines. The heating contours show greatly improved spanwise 
uniformity in Fig. 7(b) in the interior. Some end effects are evident at the top and bottom boundaries. The 
stagnation region heating is in good agreement with the benchmark but the distribution is flatter and the 
magnitude exceeds the benchmark by approximately 25% at the peak shear location. The shear contours 
in Fig. 7(d) show good spanwise uniformity as well. However, peak shear exceeds the benchmark value 
by approximately 20%. If nodes near the span boundaries are eliminated from consideration, the spanwise 
variation is nearly zero in every surface quantity as shown in Fig. 7(c). 

The cause of differences between the Option 1 algorithm and the benchmark are unknown at this time. 
Work is ongoing to resolve this difference and additional comments regarding this disparity will be offered 
in the Summary of Numerical Results section. 



(a) Values of heating - blue circles, shear - red triangles, and (b) Heating contours on surface showing significant improve- 

pressure - black squares at nodes on surface. Figure shows ment in spanwise symmetry. Some edge effects at spanwise 

results for all 11 surface nodes at a given x location. boundaries still persist. 



(c) Values of heating - blue circles, shear - red triangles, and (d) Shear contours on surface showing interior spanwise sym- 

pressure - black squares at nodes on surface. Figure shows metry. 

results for 7 interior surface nodes at a given x location to 
remove edge effects. 

Figure 7. Simulation results for 7 = 7/5 and MW = 28.8 using Option 1 reconstruction on the challenge grid 
with ten spanwise cells. 
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4- Option 2 Results on Biased, Tetrahedral Grid 

Option 2 results on the test problems do not show any significant difference in solution quality. Even though 
each individual flux component was constructed from nodal values using the STVD formulation (without 
any intermediate averaging steps on primitive or conserved variables to virtual nodes as in Option 1) the 
net integrated flux around the surface of the control volume would still admit undershoots in temperature 
across the captured bow shock. The operation count for Option 2 is more expensive than for Option 1; 
consequently, without any improvement in solution quality Option 2 is not developed any further herein. 

C. Challenge Problem - Biased, Tetrahedral Grid on 3D Sphere 

The sphere test problem 1 uses Poo = 4167 m/s, = 0.0216 kg/m 3 , = 300 K, and T wa u = 800K. 
Three sets of figures are presented to illustrate solution quality using the baseline, quasi-one-dimensional 
reconstruction on hexahedral and tetrahedral grids and using the multi-dimensional reconstruction (Option 
1) on tetrahedral grids. The three sets include surface pressure in Fig. 8, surface heating in Fig. 9, and 
surface shear in Fig. 10. A close-up view of the hexahedral surface grid and tetrahedral surface grid are 
presented in Fig. 11. 

The three sub-figures on the left of Fig. 8 from top to bottom show that surface pressure contours are 
consistently predicted on both hexahedral and tetrahedral grid systems with either reconstruction algorithm. 
The axisymmetry is well maintained with all contour lines appearing as concentric circles. On the right side of 
the figure an axisymmetric solution generated with the structured grid code LAURA is used as a benchmark 
shown as a solid black line. The symbols indicate the 3D unstructured solution cuts at 0 degrees (along the 
x-axis) and at 45 degrees so that slices cut the surface grid at two different angles. The cut lines repeat 
across the axis. In the case of perfect symmetry, all symbols lie on top of each other. In the case of imperfect 
symmetry two sets of blue squares and two sets of red squares may be observed at any radial location. The 
three sub- figures on the right indicate good agreement with the structured grid benchmark and nearly perfect 
overplotting indicating excellent symmetry. Surface pressure is well predicted with every permutation of grid 
and reconstruction algorithm tested here. This result is consistent with results obtained on the spanwise 
cylinder test. 

The solution quality for heating and shear are quite different than for pressure. The three sub-figures on 
the left of Fig. 9 show that the predicted heating in the stagnation region is sensitive to local grid orientation 
even in the case of a hexahedral grid system (Fig. 9(a)). A cross shaped contour pattern is formed in the 
stagnation region which transitions to concentric circles at approximately 40% of the maximum radius. 
(Fig. 9(b)) The standard reconstruction on the tetrahedral grids produces an even worse heating pattern as 
judged by lack of symmetry (Fig. 9(c)) and comparison to the benchmark solution in Fig. 9(d). These results 
are consistent with previous investigations. The Option 1 results with 3D reconstruction recover significant 
improvement in symmetry. A jag in the heating contours persists on the left side of Fig. 9(e) associated with 
a discontinuous change in the formation of the diagonal from the original hexahedral grid. Comparison with 
the benchmark solution at the stagnation point is improved compared to the ID reconstruction on both the 
hexahedral and tetrahedral grids (Fig. 9(f)). However, the heating level drops below the benchmark beyond 
the 60% location of maximum radius. Note that the hexahedral grid system retains excellent comparison 
with the benchmark across this same range (Fig. 9(b)). 

The quality of the predicted shear follows much the same pattern as the quality of the predicted heating 
as measured by symmetry (concentric circular contours) and comparisons to the benchmark. The hexahedral 
results show good symmetry and compare well with the benchmark (Figs. 10(a) and 10(b)); however, shear 
goes to zero at the stagnation point so it is not a sensitive indicator of solution quality in the stagnation 
region. The ID reconstruction on tetrahedral grids provides a very poor simulation of shear by metrics of 
symmetry and difference from the benchmark. (Figs. 10(c) - 10(d)). The 3D reconstruction shows good 
symmetry but overpredicts the peak shear by approximately 8% relative to the benchmark. The shear peaks 
later than the benchmark and remains larger than the benchmark values beyond 70% of the maximum radius. 

Two cuts are made through the stagnation streamline of the shock layer at 0 degrees and 45 degrees to 
show temperature contours (Fig. 12). Temperature contours from the ID reconstruction (Fig. 12(a)) do not 
overlay in the stagnation region but otherwise show expected symmetry. The temperature contours from 
the 3D reconstruction (Fig. 12(b)) show excellent symmetry even in the stagnation region. 
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(a) Hex Grid, Surface Contours 


(b) Hex Grid, Line Cuts at 0 and 45 degrees 
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(c) Tet Grid, ID reconstruction, Surface Contours 


(d) Tet Grid, ID reconstruction, Line Cuts at 0 and 45 de- 
grees 
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(e) Tet Grid, 3D reconstruction, Surface Contours 


(f) Tet Grid, 3D reconstruction, Line Cuts at 0 and 45 de- 
grees 


Figure 8. 


Surface pressure on sphere test problem. 
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(c) Tet Grid, ID reconstruction, Surface Contours (d) Tet Grid, ID reconstruction, Line Cuts at 0 and 45 de- 

grees 




(e) Tet Grid, 3D reconstruction, Surface Contours (f) Tet Grid, 3D reconstruction, Line Cuts at 0 and 45 de- 

grees 


Figure 9. Surface heating on sphere test problem. 
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(c) Tet Grid, ID reconstruction, Surface Contours (d) Tet Grid, ID reconstruction, Line Cuts at 0 and 45 de- 

grees 



(e) Tet Grid, 3D reconstruction, Surface Contours 



(f) Tet Grid, 3D reconstruction, Line Cuts at 0 and 45 de- 
grees 


Figure 10. Surface shear on sphere test problem. 
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(a) Hex Grid 



(b) Tet Grid 


Figure 11. Surface grids on sphere test problem focused near stagnation point. 
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(a) Tet Grid, ID reconstruction, Shock Layer Con- 
tours 


(b) Tet Grid, 3D reconstruction, Shock Layer Con- 
tours 


Figure 12. Shock layer temperatures on slices at 0 and 45 degrees on sphere test problem. 


D. Summary of Numerical Results 

The option 1 results presented here requires three reconstructions per element whereas the baseline ID 
algorithm requires one reconstruction per edge. The cylinder challenge problem with 10 spanwise cells 
required 14.78 seconds per relaxation step as compared to 6.84 seconds per step for the baseline ID algorithm. 
Some improvement is expected when the inviscid and viscous terms are computed in a single loop over 
elements. This refactoring will eliminate some redundancies occurring in both the inviscid and viscous 
formulations. 
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The baseline ID algorithm provides good simulation of pressure but poor simulation of heating and shear, 
with as much as ±20% dispersion about a symmetric distribution. The 3D reconstruction algorithm provides 
good simulation of pressure and eliminates almost all dispersion about a symmetric distribution on these 
biased grid test problems. The computed heating rates and shear stresses still show some differences with 
respect to the benchmark structured grid solutions. The comparisons with the sphere are better than the 
comparisons with the cylinder. The stagnation region heating for the sphere is within 4% of the benchmark. 
The stagnation region heating for the cylinder is within 8% of the benchmark. The peak shear on the sphere 
is within 8% of the benchmark. The peak shear on the cylinder is within 20% of the benchmark. There 
is no clear explanation why the sphere simulation is better than the cylinder simulation. The free stream 
conditions are different and the normal distribution of the original structured grids are different, with fewer 
nodes and greater stretching in the cylinder test problem. Certainly additional grid convergence tests are 
required here but the central point of this comparison is to evaluate the quality of an unstructured grid 
simulation relative to an equivalent (identical nodal distribution) structured grid distribution. In this regard 
it is evident that the 3D reconstruction (option 1) retains symmetry where expected in heating and shear on 
biased grids but does not necessarily attain a grid converged simulation with equivalent nodal distributions 
of a hexahedral grid system. 


V. Concluding Remarks 

The ultimate goal for any flow field simulation is to achieve a grid converged answer to a user specified 
target function (lift, drag, heat load) within a user specified uncertainty in a fully automated algorithm. Au- 
tomated grid adaptation through localized enrichment, coarsening, and movement is a necessary component 
of such a simulation capability. The present work is motivated by a vision that an underlying unstructured 
grid system composed of simplex tetrahedra in three dimensions provides the greatest flexibility to adapt 
to shocks and shear layers in a hypersonic flow. However, it is recognized that hypersonic simulations on 
such a grid system produce poor quality heating and shear - associated with algorithms that are sensitive to 
element edges skewed to critical flow structures. The present paper introduces a three-dimensional flux re- 
construction algorithm that is less dependent on edge orientation and produces improved results for heating 
and shear, even on highly biased, tetrahedral grid systems. 

Two test problems have been investigated: a hypersonic flow over a three-dimensional cylinder, including 
resolution in the spanwise direction and hypersonic flow over a three-dimensional sphere. The tetrahedral 
cells used in the simulation are derived from a structured grid where cell faces are bisected across the 
diagonal resulting in a consistent pattern of diagonals running in a biased direction across the domain. This 
grid is known to accentuate problems in both shock capturing and stagnation region heating encountered 
with conventional, quasi-one-dimensional inviscid flux reconstruction algorithms. The test problem therefore 
provides a sensitive response to algorithm on heating. 

The three-dimensional reconstruction shows significant improvement in symmetry (less than 3% disper- 
sion) of heating with some penalty of ringing at the captured shock. While it retains symmetry where 
expected in heating and shear on biased grids it does not necessarily attain a grid converged simulation with 
the equivalent nodal distribution of a hexahedral grid system. Cost for inviscid reconstruction is almost a 
factor 3 more work than the baseline algorithm using equivalent storage requirements. The net penalty on 
the complete simulation is just over a factor of two in the cylinder test problem. Some cost can be offset by 
the ability to combine the inviscid and viscous formulations in the same loop to make best use of data in 
cache. 
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VI. Appendix A - Three-Dimensional Algorithm 


The flux components are defined using nodal weights based on the Green-Gauss formulation for derivatives 
in the component directions. Thus, 


di a 
dx' 


* ’a ' 


k— faces 


(28) 


where ffc is the average value of f across all nodes defining surface k of element A. Following in parallel the 
development of the two-dimensional algorithm and defining node 4 as the additional node required to define 
the tetrahedral element one obtains 


— ^3,x'(^2,x' + fl,*' + ^4,x') + CH,x' (^3,x' + ^2,x' + f 4,x ') 

+ O 2,x'(fl,x' + ^3,x' + f^x') + <^4,x'(fl,x' + t 2,x ' + f3,x') 


(29) 


where 


Ofc,x' 


Akrik,x' 

3Q A 


and k is a face index identified by the number of the opposite node. The sum Ylk=faces a k,x' 
closed element. It is convenient to reorder and scale the coefficients as follows. 


(30) 


0 for any 


-foJ — ( 0-2, x ' + a 3,x’ + a 4,x')fl + ( a 3,x' + a l,x' + a 4,x’)^2 
+ (oi,x' + a 2, x' + &4 :X >)f3 + (oi^' + 02 , 2 ' + 03 , x ')f 4 

Fa = ~Ax' '^ 2 > a: '^ 2 + 03,x'^3 + P4,x’f-4\ 

Note that Ax' can be interpreted as a distance across the element in x' direction. It is defined 


1 

Ax’ 


V! (\ a n-l,x' 
n=nodes 


^n+ l,x' “1“ ^n+2,x / |) 


(31) 

(32) 

(33) 


where a cyclic indexing is assumed. The reordering and scaling yield the following relations for /3„ )X /. 



V! I Pn,x' I 


Ax (o n _ i,^;' -(- cy nA i x f o n _|_2,x') 

0 

1 


(34) 

(35) 

(36) 
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A. Option 1 - Virtual Node Averaging in 3D 

The equations developed in the original section for two-dimensions are unchanged for three dimensions except 
that the sum over nodes includes a fourth node required to define a tetrahedral element. 


B. Option 2 - Weighted Average of Edges to Principal Node in 3D 

Select the largest of \0 n ,x'\ to identify the principal node for construction of f x j and express it as a function 
of the remaining coefficients. Assume node 3 is the principal node. Then 


dfx' = 01,x'fl + 02,x'h + 03,x'f3 + 04,x'f4 

= 01,x'fl + 02,x'f2 + 04,x'f4 — (01, x' + 02, x' + Pi,x')h ( 37 ) 

= 01, x' ( fl ~ f 3 ) + 02,x'{f'2 — /3) + 04, x' (/*4 — f3) 

The reconstructed flux in the direction x', f x >, is computed as a weighted average of surrounding edges. 
Furthermore, it is noted that if any of the surrounding edges are parallel to x' then the weight of that edge 
will equal 1 and the weight of the other edges will equal zero. 


fx‘ — \01,x‘\^l-3,x' + \02,x'\^2-3,x' + \04,x'\^4-3,x' 


(38) 


where 


fl-3,x' 
^2-3,x' 
f.4— 3,x' 


1 - 

2 . 
1 - 

2 . 
1 - 
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fl,x 

f 2 ,x 


*4,x 


' + f 3 , 
' + f 3 , 
' + f 3 , 


x' 


x' 


x' 


sign (/?i ,x' ) R]"- 3,®' I- A-i- 3,®' I (dqi-3,*' 

sign(/?2,x') R 2"-3 ia: ' I A 2-3 ,x' I ( dq_2-3,x ' 
sign (/? 4 , s' ) R 4- 3 ,x' I A 4-3,x' I (rfq4-3 >a; ' 


dc\\—3,x',lirn) 

dd\2—3,x',lini) 

dd^4—3,x',lini) 


dq 1-3, x' 

dq 2-3, x' 
dq_4-3,x' 


dqi — 3,x' ,lim 
dq2—3,x , ,lim 
dci4-3,x' ,lim 

d(\n,x' 


Ri-3,x'(qi - q3) 
R2-3,x'(q 2 — q3) 
R 4 -3,x'(q4 - q3) 


minmod 


2dqi,x', 2dqi_ 3 ,x', 2dq 3>a: /, 


minmod 


2dq 2 ,x' , 2dq 2 _3,x' , 2dq 3 ,x' , 


minmod 


2dq4,x' 5 2dq4_3 >a: ' , ^dq^^' , 


^(dqi, x ' + c?q3,x') 

^(dq- 2 ,x' + dq 3 ,x’) 
^(dq4,x' +dq 3 ,x') 


Ax R/; ,,f ‘ Vf J _ / .3,p 1,,'' 


(39) 

(40) 

(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 


C. Component Assembly on Faces in 3D 

The flux across each equivalent face An in 3D is constructed from the component flux in the principal 
direction and orthogonal to the principal direction. 

^ An = HAn,x'^x' T ^An,y'^y' T ^An,z'^z' (49) 

where An refers to face separating node n from other nodes in the element. 
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